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We solve numerically the Boltzmann equation describing the evolution of a cosmic string network 
which contains only loops. In Minkowski space time the equilibrium solution predicted by statistical 
mechanics is recovered, and we prove that this solution is stable to non-linear perturbations provided 
that their energy does not exceed the critical energy for the Hagedorn transition. In expanding 
Einstein - de Sitter Universes we probe the distribution of loops with length much smaller than 
the horizon. For these loops we discover stable scaling solutions both in the radiation and matter 
dominated epochs. The shape of these solutions is very different in the two eras, with much higher 
energy density in the radiation epoch, and a larger average loop length in the matter epoch. These 
results suggest that if the conditions for formation of loop networks are indeed satisfied, these could 
in principle be good candidates for structure formation. 
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I. INTRODUCTION 

One of the greatest challenges of modern cosmology is 
to explain the large scale structure of the universe: how 
did the universe, which is known to have been very ho- 
mogeneous at early times, eventually develop structures 
such as galaxies and anisotropics in the cosmic microwave 
background (CMB)? There are currently two different 
classes of theory which try to provide an explanation — 
inflation and topological defects g| — though the two 
may also be combined (^-^]. In inflationary scenarios, ini- 
tial vacuum quantum fluctuations eventually lead to the 
formation of structures. In topological defect theories a 
network of vacuum defects is produced at a phase tran- 
sition in the early universe. These defects naturally pro- 
vide a source of inhomogeneity for all subsequent times. 

Although there are several classes of defects , cur- 
rently the most popular candidate for structure forma- 
tion is cosmic strings. These are line-like vacuum defects 
fully specified by a single parameter, fi, their energy per 
unit length ||-|l3| . For strings to be viable candidates for 
structure formation, they should neither come to domi- 
nate the energy density of the universe nor disappear al- 
together. Instead they should "scale" , by which we mean 
that the energy density in strings should remain a fixed 
fraction of the energy density in the universe. 

The formation and evolution of a network of cosmic 
strings is extremely complex. A large amount of analyt- 
ical ]14|-[l7],^H and numerical ]l9|- pl|Jl8t work has been 
devoted to exploring the evolution of these networks in 
expanding flat universes with no cosmological constant. 
Current wisdom is that the network quickly evolves to- 
wards a fixed scaling solution, though the details of this 
solution differ between the radiation and matter eras. In 
this attractor solution the network is dominated by infi- 
nite strings. 

In a widely known but unpublished manuscript, Fer- 
reira and Turok p3{ were the first to conjecture that a 
different scaling solution might be reached if only string 
loops were formed during the phase transition. In this 



different scaling solution no infinite strings are present, 
given that none exist in the initial conditions: the net- 
work is forever made up solely of loops. This conjecture 
is an intriguing possibility, with potential far reaching 
consequences for structure formation. Undoubtedly the 
time-time correlators of the stress energy tensor of such 
networks is rather different. Radically different predic- 
tions for galaxy clustering and CMB anisotropics could 
ensue. 

It was later argued that networks containing only loops 
might indeed form after first order string producing phase 
transitions |24| - |26f| . Such transitions proceed via bubble 
nucleation, and as the bubbles expand and collide, string 
loops are produced. It was suggested that this mecha- 
nism may be vastly more effective than the more conven- 
tional Kibble mechanism. 

In this paper we use the model of Ref. to study 
the evolution of networks of string loops through a 
Boltzmann-type rate equation. Previously this equation 
was studied analytically in detail in references p2| , p7j] . 
However, given the complicated nature of the equation 
it is a non-linear and non-local, integro-partial differen- 
tial equation — only asymptotic solutions were searched 
for, that is, the focus was on the behaviour of loops 
whose length was much bigger than the horizon size. It 
was argued |^2| that a stable scaling solution exists in 
Minkowski space time as long as the initial energy den- 
sity in loops is less than a critical value which corresponds 
to the Hagedorn energy. In a radiation dominated uni- 
verse, a scaling solution was also found. However, this 
was thought to be unstable — depending on the initial 
energy density in loops, these were either thought to dis- 
appear, or infinite strings were predicted to form. In the 
matter era, no stable solution was found, and for much 
of parameter space the loops were expected to disappear. 

The purpose of this paper is to search for numerical 
solutions of the full rate equation of p^j, rather than 
asymptotic solutions. This is a challenging numerical 
task given the complicated nature of the rate equation, 
and in fact we are only able to probe accurately the be- 
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haviour of small loops. The reasons for this and some 
of our methods for tackling the numerical problems are 
described in the rest of this paper. We also describe the 
results of our code and elaborate on the properties of the 
solutions found. 

This paper is organised in the following way. In Sec- 
tion [n] we review the rate equation and set up our no- 
tation. In Section III we discuss some further analytical 
massaging connected with the practical aspects of our 
numerical implementation. The core of our paper is in 
Sections IV and |v|, where we present our numerical solu- 
tions, as well as some analytical results used as a check 
upon the numerics. The results presented there can be 
summarised as follows. In a Minkowski background we 
recover the scaling solution predicted by statistical me- 
chanics. We prove that this solution is stable against non 
linear perturbations, as long as these do not exceed the 
critical mass for the Hagedorn transition. In expanding 
Universes and for loops much smaller than the horizon, 
we discover stable scaling solution both in the radiation 
and matter dominated epochs. The shape of the solution 
is quite different in the two eras, with much higher en- 
ergy density in the radiation epoch, and a larger average 
loop length in the matter epoch. Based on these results, 
we suggest that if the conditions for formation of loop 
networks are indeed satisfied, these could in principle be 
good candidates for structure formation. Finally, in Sec- 
tion VI we conclude, outlining further work prompted by 



our results. 



II. THE RATE EQUATION 

In this section we set up our notation and review 
the rate equation proposed in for the evolution of 
networks of loops. For clarity we introduce the non- 
expanding and expanding universe cases separately. The 
rest of the paper will focus on searching numerically for 
solutions of the equation. 



A. Non-expanding Universe 

Consider first the evolution of a network of cosmic 
string loops in a non-expanding universe. In this case one 
would expect the network to reach an equilibrium, time- 
independent distribution. We will study the approach to 
equilibrium and the propties of that solution. 

Let n(£, t) be the number of loops per unit volume with 
physical length between £ and £ + d£ at a time t. In p^ j 
it was argued that in a non-expanding Universe n(£, t) 
satisfies 
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A{1] £', £ - £')n(£', t)n{£ - £', t)d£' 
B{£JJ -£;£')n(£',t)d£' 



n(£,t) / A(£ + £';£,£')n(£',t)d£' 
Jo 

f £/2 

n(£,t) / B(£',£-£'-£)d£', 
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(1) 



where the four integrals describe collision interactions be- 
tween the loops. The coefficient A(£ + £';£,£') gives the 
probability that two separate loops of length £ and £' 
collide and produce a loop of length £ + £' (assuming as 
usual that in the collision partners are exchanged), and 
it is given by [^2| 



A(£ 1 +£ 2 ;£i 1 £ 2 ) = xh£2, 



(2) 



where x ' is a relative velocity. The function B (£',£ — £'; £) 
gives the probability that a loop of length £ self-intersects 
to produce two loops of length £' and £ — £' . Arguments 
in p2|j27|] suggested that 



B(y,£-y;£) = 



X^-y)(£ + Q 5/2 

e' 2 {£-v + C) 5/2 (y + C) 5/2 ' 



(3) 



where % is another relative velocity. The two lengths in 
(||) are £, the distance along which the Brownian loops 
are correlated in direction, and the distance between 
any kinks on the loops. These lengths are assumed to be 
time independent in a non-expanding universe. 
The total length density of loops, L is given by 



(4) 



L(t) = / n(£,t)£d£ 



and in Minkowski spacetime this remains constant since 
energy is conserved. In fact the proof of this, substiting 
the rate equation into the time derivative of (|J) is quite 
complex as one often has to deal with diverging integrals 
p7[ . This problem of infinities will also have to be dealt 
with in our numerical work, where we shall use (^) as a 
check on our code. 

The total loop number density is given by 



N(t) 



n(£,t) d£. 



(5) 



There are no limitations on this quantity, but combined 
with the energy E ~ L it gives information on the shape 
of the loop distribution. 



B. Rate equation in Expanding-Universe 



In an expanding Universe the rate equation becomes 



dn ,„ ,d£ dn(t,t)- ,. 
at o£ a£ e^o+ 



n{i,t)£ 8{£) 



- 3Hn(£,t) 
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+ / A(£;£',£-£')n(£',t)n(£-£',t)d£' 
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/ B(£,£' -£;£')n{l',t)dl' 
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n{£,t) / A(£ + £';£,£')n(£',t)d£' 



n{£,t) / B(l',e-£';£)dl', 



(6) 



where t is the rate of change of the length of a loop. This 
has contributions from expansion, redshift and gravita- 
tional radiation, and may be written as p2] 



Kl 

T' 



(7) 



where we expect < K<0.1 [|22|]27]]. The third term 
in @ guarantees that n{(.,t) — 0, W < Vi, whilst the 
fourth term is the contribution from the expansion of the 
universe with H being the Hubble parameter. The func- 
tions A and B in (Q) are the same as those of equations (|^) 
and (||) respectively, where we again assume the same two 
length scales on the loops, £ (distance along which loops 
are correlated in direction — c.f. 3-scale model |TJ]) and 
£ (the distance between kinks — c.f. fll4||). Now though 
these two scales are assumed to 'scale' 



(8) 



where £ c and £ c are both constants. As discussed in |2S 
and in analogy with p4| we expect 



(9) 



Our aim will be to see whether the assumption (||) 
combined with equation (||) are consistent with the en- 
ergy density in loops also scaling. For that reason it is 
convenient to rewrite the rate equation (^|) in the scaling 
variables 



_ £ 
Z ~ t 
g(z 7 t) = t i n{e,t) 7 



(10) 
(11) 



so that scaling energy density corresponds to g(z, t) = 
g(z). Recall that from the definitions of A (^) and B (||) 
one has 

A(zt; z't, (z - z')t) = t 2 A(z; z', z - z'), (12) 
B(zt; {z — z)t, z't) = —B(z, z — z; z ) (13) 

where it is important to remember that in B(z, z' — z; z'), 
£ must be replaced by £ c . Hence the rate equation for 
g(z, t) becomes: 



dg{z,t) 



g{z, t) [ 4 
dg(z,t) 



3 --K 
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dz 



(z(1-K)+TGij,) 



TG^(z)[g(z,t)]\ z _ 0+ 
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+ I A(z;z',z - z')g(z',t)g(z - z',t)dz' 



+ / B(z,z' - z;z')g(z',t)dz' 



g(z,t) A(z + z':z,z')g(z',t)dz' 



j.z/2 

g(z,t) / B(z! ', z — z' '; z)dz' '. 
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(14) 



Here the parameter p is defined by the expansion rate; 
R oc t 1 ' p , so that p = 2 in the radiation era, and p = | 
in the matter era. Observe that the RHS of the above 
equation is not explicitly t dependent and hence that it is 
consistent to have g(z,t) = g(z) - i.e. a scaling solution 
is consistent with the equation. 

Our aim is to solve equation (|l4|) to check whether such 
a scaling solution is indeed reached, and if so what are 
its properties. 



III. NUMERICAL IMPLEMENTATION 

A. Variables suitable for numerical work 

Any numerical approach to solve equation ( |l4| ) has to 
first deal with the obvious numerical problem posed by 
the semi-infinite interval [0, oo) in which the z variable 
lies. As was discussed in p2| , p7| , just chopping off the in- 
tegrals in jT4] ) at some upper value of z is not a resolution 
to this problem as this changes completely the nature of 
the solutions of the rate equation. Instead we introduce 
a new variable: 



1 + z/c n 



(15) 



in which c max is a fixed number which will be chosen to 
optimise our numerical analysis. Note that, as opposed to 
z, the variable y now lies in the finite range < y < c max . 
Also it will be useful to write down the inverse of ([151); 



z (y) = 



u 



1 - y/cn 



We then define 



Hy( z ),t) = g( z ,t) 



(16) 



(17) 



so that the boundary condition g(z — » 00, t) — Vi be- 
comes 



h(c max ,t) = Vi, 



(18) 



and may now be imposed without loss of generality. This 
boundary condition is necessary so that the energy den- 
sity in loops is finite. 

Furthermore, in an expanding universe only, we change 
variables to logarithmic time 
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t — > u — ln(i), 



(19) 



so that on letting h(y,u) = h(y,t(u)) the rate equation 
becomes 



<%(y, u) 



(1 - y/c max f (z(l - K) + TGfi) 



dy 

, Uj\ |y — >0 + 



+ / - A(z(y'),z(y)-z(y'))h(y',u) 
Jo 

dy' 



xh(y(z(y) - z(y')),u) 



(1 y / Cffiax)^ 

I h(y',u)[B{z(y),z{y'))-h(y,u) 
dy' 



x A(z(y),z(y'))] 



(1 ~ y'/c m ax) 2 



hi,, a ) f A(z(y), %'))%', u) — £f ^ 



-h(y,u) / B(z(!/V(!/))7] 7 



dy' 



(20) 



Note that for clarity we have now omitted the unneces- 
sary arguments of the A and B coefficients: for exam- 
ple, we have written A(z; z' , z — z') = A(z' , z — z') and 
B{z,z' - z;z') = B(z,z' - z). 

In a non-expanding universe the rate equation is that 
given in (^0|) but where one retains only the four integral 
terms on the RHS, and where there is no logarithmic time 
(so t = u). 



B. Choice of parameters 

In order to solve ( |20| ) numerically, we clearly have to 
choose values for the various parameters which appear in 
the rate equation as well as a value for c max . 

Firstly, by comparison with numerical simulations of 
cosmic string evolution in expanding universes (see the 
discussion in section 3 of |E2[ and also E3 ) , we have taken 



X = X = 0.2 



(21) 



though in general these relative velocities need not be 
equal. The values given in ( |2l| ) were also used for our 
simulations of loops in Minkowski space. 

The correlation length £ c of equation (|J) was taken to 

be 



& = 1 



(22) 



which is of the order of the horizon size. Note that the 
only place in which £ c appears in the rate equation is in 
the coefficient B[z' , z; z + z'), where the relevant combi- 
nation on constants is x/^ 2 ■ Thus we can change the 



size of £ c simply by scaling \- The value of £ c is approx- 
imately given by (|9j), so we have taken 



C c = 5.0 x 10" 



(23) 



since for GUT cosmic strings G\i ~ 10~ 6 , and 50<r<100 
(we have chosen T = 50). In Minkowski space, we took 
£ = 1 and ( — 5.0 x 10~ 5 consistenly with the assump- 
tion of the time- independence of these variables in a non- 
expanding universe. 

As mentioned above and discussed in [22) , < K < 0.1 
and we have taken 



K = 0.05. 



(24) 



We will also consider the effect of varying K about this 
value, since it was shown in analytical analysis of |2^ j 
that the precise value of K may influence whether or not 
a scaling solution can exist in the matter era. 

Finally, observe that the probability B(z' , z — z'; z) is 
very strongly peaked about z' = z — z' £ c . In order to 
deal with such rapidly changing functions in the region 
of z = C c we found that greatest numerical accuracy was 
achieved by choosing 



10" 



(25) 



We then divided the < y < 10~ 4 interval into N equally 
spaced intervals (with N = 200 or N = 500), and used 
open integration routines and adaptive mesh techniques 
in the u direction to solve equation (20). 

We end this section with the following important ob- 
servation. By choosing c max = 10 , our numerical work 
can deal very accurately with smallloops, but cannot deal 
accurately with large loops. The reason is the following. 
Let us assume that we work with N = 500 intervals in 
the y direction. Then the last interval corresponds to 
9.998 x 10~ 5 < y < 10 -4 , or alternatively using equation 
( |l6| ) to 0.0499 < z < oo. In other words, the last inter- 
val in y contains all the loops with length z>£ c /50, and 
the evolution of these loops in time can clearly not be 
followed accurately using only one interval. One might 
then be concerned that inaccuracies with loops of length 
z>£ c /50 might lead to inaccuracies in the evolution of 
those with smaller lengths since the rate equation (|l4|) is 
non-linear as well as non-local. This indeed will be true 
though it is difficult to estimate the size of these effects. 
However, an expansion to lowest order in z/£ c of the rate 
equation ( |l4"| ) shows that the small z behaviour decouples 
from the large z one ETj (though this is no longer true 
to higher order in z/Q^T 

We therefore conclude this section by saying that our 
numerics will be probing the small z behaviour of the loop 
distribution g(z, t) and not the large z ^> Cc Cc one which 
was the subject of p2]. This numerical study is therefore 
complementary to the analytic one of |p2| rather than 
being at odds with it. 
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IV. LOOPS IN MINKOWSKI SPACE TIME 

Given these comments, we begin by considering the 
rate equation ([!]) in a non-expanding universe. (As com- 
mented above for the purpose of the numerics we work 
with ( p0| ) where only the four integrals on the RHS are 
kept and where u = t). 



A. Some analytical results 

In Minkowski space we know the equilibrium solution 
of (Q): combining the principle of detailed balance with 
the form of the A and B coefficients gives the equilibrium 
(and hence stable) distribution as |22j 



:ime=1 ,2000, 4000,..., 1 10000 



n(£) 



-01 



-01 



x £3/2 {l + C) 5/2 



(26) 



Substitution of ( pq ) into (|I|) shows that this is indeed a 
solution with (3 arbitrary. However, the energy density 
in loops must be finite so that /3 > 0. 

The solution (^6|) is consistent with the statistical me- 
chanics prediction that in equilibrium [E8| 



i{£) = ce 



9-5/2 



(27) 



for £ £. It has also been shown [£8| that at (3 — when 
the loops have a scale invariant distrubition, the energy 
density in loops is equal to a critical value and a phase 
transition occurs. This is the Hagedorn transition for 
which the loops network is unstable against the formation 
of an infinite string. 

In [P2l , the approach to the equilibrium solution (|26| ) 
was studied using the rate equation ([!]). This was done 
by substituting distributions of the form 



n(£,t) 



-0t 



Cn(t) 



•C) n 



(28) 



with c„ and (3 are arbitrary parameters, into the rate 
equation. Equations for c„ and (3 were then obtained by 
studying the £ ^> C limit of the rate equation. Analysis of 
these equations showed that a for an initial distribution 
of the form (p8|), a stable equilibrium solution is reached 
which is characterised by 



Co 



Cn>l 



0. 



(29) 



Furthermore, since energy is conserved by the rate equa- 
tion, (3 was shown to change as well, with its equilibrium 
value such that the energy density at equilibrium is the 
same as the initial energy density. This evolution was 
shown in figure 1 of fl22|| . 

If the initial energy density in loops is larger than the 
critical energy density, then one would expect from the 
results of statisical mechanics that infinite strings are 




FIG. 1. Initial conditions with co = 0.8 < c eq , c n >o = 0, 
f3o = 0.01. We have plotted the function n(£,t) for times from 
t = 1 to t = 110000. 



formed. A plausibilty analysis for this (based on the ap- 
pearance of non-analytic terms when (3 — > 0) was given 
in |22] | . In the present paper we are concerned only with 
distributions that have a total energy density less than 
that for the Hagedorn transition. Analysis of the case 
with energy greater than the critical one will be given in 
a future publication p9f (see also |30|). 



B. Numerical results 

We have solved ([!]) numerically, starting with a vari- 
ety of initial conditions parametrised as in (|2^) , but not 
necessarily close to the equilibrium configuration, though 
always with subcritical energy density. 

Figures ([!]) and (||) shows the time evolution given the 
initial conditions, cq — 0.8 and cq = 1.2, respectively, 
with (3 = 0.01 and c„>o = in both cases. Recall that 
with our choice of parameters c eq = 1. In the first case 
Co approaches c eq from below, and in the second one c eq 
is approached from above, whilst in both cases the power 
dominant power law behaviour is indeed (I + C)~ 5 ^ 2 as 
expected from ([26]). (In other words, all the c„>o remain 
zero for all times.) It is difficult to see from these plots, 
however, if (3 changes as predicted since for these small 
values of £, e' 131 = 1 to very high accuracy (see further 
comments on this below). 

The only point of the distribution which does reach 
the scaling solution is n(0,i), which in fact does not 
change. It can be shown analytically from (|l]) that the 
time derivative of n at £ = is indeed zero. There is very 
little meaning in loops of size zero, so we do not comment 
further on this exception to the equilibrium distribution. 
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time=1 ,2000, 4000,. ..,150000 




In figure the initial conditions have cq = 0, ci = 1.2, 
Cn>i = and [3 — 0.01. That is, this initial distribution 
has a power law behaviour of (£ + C)~ 7 ^ 2 - Despite that, 
the plot shows that as the network evolves, c\ — > and 
Co — * 1 so that the system evolves to its equilibrium dis- 
tribution (though again it is difficult to comment on (3). 

The stability of this solution with Co = c eq and c„>o = 
is very important. Previous proofs of stability have 
been confined to linear perturbations. No analytical evi- 
dence has been given for stability in the nonlinear regime. 
Our results bear strong evidence that the solution is 
reached with any initial conditions, and thus that it is 
stable also in the non-perturbative regime. 

An important test for the code is to check whether the 
total loop length density (energy) stays constant. The 
simulations show that the more sampling points (the 
larger N), the better the conservation law is satisfied. 
We decided to accept sampling densities leading to mass 
variations of ~ 5%. 



V. EXPANDING UNIVERSES 



FIG. 2. Initial conditions with Co = 1.2 > c eq , c n >o = 0, 
P = 0.01. We have plotted the function n(£, t) for times from 
t = 1 to t = 110000. 



time=1 ,2000,4000,...,! 300000 




FIG. 3. Initial conditions with Co = 0, Ci = 1.2 > c eq , 
c n >i = 0, P = 0.01. We have plotted the function n(£,t) for 
times t = 1 to t = 130000. 



A. Analytical results 



In 1 22 1 the asymptotic behaviour of solutions to the 
rate equation (jlj) were searched for. More specifically, 
using the ansatz 



g(z,t) = 6(z) 



(z + 



Cn{t) 



(z + cy 



(30) 



substituting into the rate equation (|14|), and assuming 
(z ^> £ c 3> Cc) led to equations for $ and c„. Analysis of 
these led to the conclusion that, at least to order n = 1, a 
scaling solution exists in the radiation era but not in the 
matter era. Furthermore, the radiation era solution was 
thought to be unstable. That is, depending on the initial 
energy density in the loops, these were either thought to 
disappear or infinite strings were thought to be formed. 
In the matter era no stable solution was found and for 
much of parameter space the loops were expected to dis- 
appear p2j . 

Before presenting our numerical analysis of this ex- 
panding universe case, we introduce the concept of the 
mean loop length since this will be useful for the analy- 
sis of our numerical solutions. Recall that by definition, 
scaling of the energy density of the loop network means 
that 



dE Si 
dt 



0, 



where 



E sc = A 4 / zg(z,t)dz. 
Jo 



(31) 



(32) 



As was shown in [E7 , £, r satisfies 
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Radiation era. Time= 1 , 1 .5, 2,. ..50 



2.0x10 r 




Time evolution 



FIG. 4. Evolution of the network in the radiation era. Ini- 
tial conditions were Co = 1.2, c n >o = and P = 10" 2 . We see 
that the distribution after 10 to 20 logarithmic times reaches 
a stable scaling solution. We have plotted the function h(y). 



dt 



f- 



(2 + K -3/p)E ac (t) -TGnN sc (t) (33) 



where N sc (t) is the total number of loops in scaling vari- 
ables, 



N sc (t)= / g(z,t)dz. 
Jo 



(34) 



As commented in |2^| , we see that the assumption of scal- 
ing lengths £ and £ does not automatically imply that the 
energy density in loops also scales in an expanding uni- 
verse: this will only happen if N sc and E sc reach values 
such that the RHS of (3^) vanishes. If that is the case, 
then the mean length (z) is given by 



(z) 



Jo°° zg{z,t)dz E sc (t) 



2 + K-3/p 



(35) 



(36) 



JZ°g(z,t)dz N sc (t) 
With our choice of parameters 



{ Z ) rad ~ i x i(r 4 = 2C C 

{z)matt ~ 1 X 1(T 3 = 20Cc 



B. Numerical Results 

1. Radiation era 

We gave the network initial conditions of the form ( |30| ) . 
but with various values for c„ and 0. In figure (|]) we 
show the time evolution of h(y, u) given the initial condi- 
tions Co = 1.2, c„>o = and — 1CP 2 . The time range 



FIG. 5. Plot of the scaling energy density, E{u) during the 
evolution of the network in the radiation era. Initial condi- 
tions were Co = 1.2, f3 = 10~ 2 and c n >o = 0. The energy den- 
sity increases rapidly just after start, and firmly, after about 
20 logarithmic times, approaches the scaling value from be- 
low. 

is u = 1,1.5, ..,50. Figure (||) shows the scaling energy 
density (defined in (|32j)) as a function of time. We see 
that the energy density after 10 - 20 logarithmic times 
approaches a stable scaling value from below. 

In figures (^|) and (|?]) we show the time evolution given 
two different initial distributions both with energies be- 
low the scaling value. It appears that these reach the 
same scaling distribution as that of figures (||) and (]?]) 
— we shall comment on the form of this distribution be- 
low. Figure (^) shows time evolution starting out with an 
initial distribution with energy density higher than the 
scaling value, Co = 12, — 0.01 and c n >o — 0. Again the 
same scaling solution seem to be approached from above. 
The same solution is also found if we set cq — but other 
of the c„>o 7^ 0. These results provide evidence that the 
small loops reach a unique stable scaling solution in the 
radiation epoch, and that this solution is an attractor 
within a very large basin of initial conditions. 

We now attempt to characterise more thoroughly the 
scaling solution. Its energy density is given by 



1.2 x 10 d 



(37) 



We could try to fit its shape with the ansatz g(z) = 
ce _/3z /(z + C) 5 ^ 2 - This is the leading order term in the 
ansatz (^) which was used in ^] for looking at large 
loops. We used a least squares method and the result 
is plotted in Figure ( ^|): clearly this function is a good 
approximation in the radiation epoch. The optimal com- 
bination of c and (3 is 



C-rad ~ 6.9, 



'rad 



34 



(38) 
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FIG. 6. Evolution of the network in the radiation era. Ini- 
tial conditions were Co = 1.2, p = 10 2 and c n >o = 0. We see 
that the distribution after 12 to 15 logarithmic times reaches 
a stable scaling solution. We have plotted the function h(y). 



Radiation era. Time= 1 , 1 .5, 2,.. .40 



•2: 2.0x10 




y 



FIG. 7. Evolution of the network in the radiation era. Ini- 
tial conditions were co = 0.2, 13 = 1(T 2 and Cn>o = 0. We see 
that the distribution after about 20 logarithmic times reaches 
a stable scaling solution. We have plotted the function h(y). 




y 



FIG. 8. Evolution of the network in the radiation era. Ini- 
tial conditions were Co = 12, (3o = 10" 2 and Cn>o = 0, which 
means we start out well above the scaling solution. We see 
that the distribution after about 20 logarithmic times reaches 
a stable scaling solution. We have plotted the function h(y). 

(Though the work of |2^] was probing a very different 
regime, there the (un-stable) scaling solution had values 
Crad = 7.1 and {3 ra d = 108.) 

Finally, the scaling solution has a mean loop length 

{z} rad, numerical ~ 1-6 X 10 (^^) 

in good agreement with the theoretically predicted value 
of (z) w 1 X 10~ 4 from equation (p6|). 

In an attempt to understand how the solution depends 
on the parameters involved, we also show results corre- 
sponding to a different value for K. Figure ( [l(]) shows 
how a scaling solution is still approached, but now with a 
lower energy density, E sc ~9x 10 2 . This is understand- 
able, since a smaller K could represent either a larger 
loss from gravitational radiation in the large loop limit 
or a greater loss from redshift of velocities. 

2. Matter epoch 

We performed a similar analysis in the matter epoch. 
The analytical analysis of p2] predicts that no scaling 
solution exists for large z in this case. Whilst this may 
well be the case, our numerical results show that a scaling 
solution seems to be reached for small z. 

Initial conditions of the form (J3TJ) were again used. 
Figure ( |l|) shows the time evolution of h(y,u) in the 
matter era given initial conditions, cq — 1.2, [3 = 0.01 
and c„>o = 0. The corresponding scaling energy density 
is plotted as a function of logarithmic time, E(u) in figure 
( |l2|). It shows how an energy density over the scaling 
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Function fitting, radiation era 

4.oxio 11 r 1 1 1 i 1 1 1 i 1 1 1 i 1 1 ^ 




FIG. 9. Attempt to fit the radiation scaling solution to 
the ansatz, g(z,t) = ce^ l3z /(z + Cf^ 2 ■ As we can see this is 
a nearly perfect fit. The constants found are, c = 6.9 and 
f3 = 34 



Radiation era. Time= 1 , 1 .5,2,... 1 00 

6x1 o 1 1 | 1 1 1 i 1 1 1 i 1 1 1 i 1 1 1 r~ 




FIG. 10. The scaling solution in the radiation era, given a 
new value for K = 0.025. We see that the solution is signif- 
icantly smaller than the previous one, with E rad w 0.9 • 10 3 
and (z) = 1.70 • 10~ 4 




FIG. 11. Evolution of the network in the matter era. Initial 
conditions were Co = 1.2, c n >o = and (3 = 10~ 2 . We see that 
the distribution immediately after the start starts decreasing 
and reshaping. The distribution goes a bit below the scaling 
solution, but catches up. The loop distribution stabilises after 
about 10-15 logarithmic times, but at value far lower than the 
radiation era. We have plotted the function h(y). 



Time evolution 
250 I — 1 — 1 — 1 — 1 — I — 1 — 1 — 1 — 1 — I — 1 — 1 — 1 — 1 — r~ 



200 - 



150 - 




100 - 



50 - 



I , , , , I , , , , I , , , , I , , , , I , , , , I 

5 10 15 20 25 

Time 

FIG. 12. Plot of the scaling energy density during the evo- 
lution of the network in the matter era. Initial conditions 
were Co = 1.2, c n >o = and f3 = 10 -2 . The energy density 
decreases rapidly just after start, but then approaches the 
scaling value from below. 
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Matter era. Time= 1 , 1 .5,2, ...30 

1.5x1 10 I 1 1 1 I 1 1 1 I 1 1 1 I 1 1 1 
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FIG. 13. Evolution of the network in the matter era. Initial 
conditions were Co = 0.1, c n >o = and /3 = 10 -2 . We see that 
the distribution immediately after the start starts decreasing 
and reshaping. The distribution then increases until it reaches 
the scaling distribution for the matter era. The loop distribu- 
tion stabilises after about 10-15 logarithmic times. We have 
again plotted the function h(y). 

value decreases rapidly, overshoots, and then increases 
to asymptote towards the scaling value after about 10 
expansion times. 

Figure ( |l3| ) shows the time evolution of h(y, u) given 
an initial distribution with an energy density lower than 
the scaling value. The function immediately starts re- 
shaping, and after some time the same scaling distribu- 
tion is found. We considered also initial configurations 
with different values c„>o ^ 0, and found that they ap- 
proach the same scaling solution. As in the radiation 
era the same scaling solution is always reached, starting 
from a large sample of possible initial conditions. Again 
though we are not very sensitive to different values of j3. 

The stable scaling solution in the matter epoch is 
rather different from the radiation epoch scaling solution. 
Its scaling energy density is nearly 10 times smaller 

£«1.5xl0 2 , (40) 

Its shape is also very different. The mean scaling loop 
length (z), is larger in the matter epoch, as predicted. 
We found 

{z} matt, numerical ~ 0.5 X 10 . (^1) 

The difference between (z) ra d and (z) m att is not as large 
as theoretically predicted, but shows the right tendency. 

Finally we tried to fit the scaling solution to an ansatz 
g(z) = ce^ l3z /{z + C) 5 ^ 2 - Again we have used a least 
squares method. The result is plotted in ( |l4| ). Clearly 
our scaling solution does not fit this ansatz. However, 



Function fitting, matter epoch 
2.0x10 I 1 1 1 I 1 1 1 I 1 1 1 I 1 1 <~ 




FIG. 14. Attempt to fit the matter scaling solution to the 
ansatz, g(z, t) = ce~ l)z / (z + () 5 ^ 2 . As we can see the functions 
are not compatible. 

it is not difficult to see that a distribution of the form 
g(z,t) = ce~ l3z /(z + £)" with n < 5/2 would fit our 
scaling distribution better. 

VI. CONCLUSIONS 

We re-examined the evolution of networks of cosmic 
strings consisting solely of loops. Our aim was to provide 
full numerical solutions to the rate equation first given in 
p2| , though in fact we were only able to probe accurately 
small loops. 

In Minkowski space time we found that for energies 
lower than the critical energy the system evolves towards 
a scaling solution. This happens for a large variety of 
initial conditions not necessarily in the vicinity of the 
scaling solution. The results obtained with our numeri- 
cal code agree with the statistical mechanics prediction, 
as well as with the analytical predictions made in p2| . 
However we also predict that the scaling solution is stable 
against non linear perturbations, as long as the critical 
mass is never exceeded. For total energies larger than the 
critical we were unable to obtain reliable results, though 
infinite strings are expected to form. 

In expanding universes and for small loops, a scal- 
ing solution was found both in the radiation and matter 
dominated epochs. The shape of the solution is quite 
different in the two eras. The gradients of the distribu- 
tions for small z agree with an analytical analysis of the 
rate equation for small z carried out in p7[ ]. The scal- 
ing energy densities are E raa - ~ 1200 and E rnatt er ~ 150, 
whereas the average scaling loop lengths, z = i/t, are 
(z rad ) « 10- 4 , and (z matt ) « 0.5 x 10" 3 . 
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Our results complement the findings in [E2| for radi- 
ation and matter dominated Universes. In both these 
cases we found evidence for scaling solutions containing 
many small loops. Thus there will always be a scaling so- 
lution in these eras, though from the results of |^] these 
may or may not contain many loops with length larger 
than the horizon. 

We close with a description of three promising research 
avenues which our work may have prompted. Firstly our 
results are very encouraging for a loop dominated sce- 
nario of structure formation. Our discovery of scaling 
solutions for small loops ensures that no cosmological 
cataclysm happens in these scenarios. Undoubtedly the 
unequal time correlators, required for computing struc- 
ture formation, are very different in loop networks and in 
the usual networks dominated by long strings. The pre- 
dicted CMB anisotropy and galaxy power spectra will 
therefore be very different. One wonders if a better fit to 
the data might be achieved in these models. 

Before this idea can be carried out, we should investi- 
gate the radiation to matter transition which is a crucial 
stage in structure formation. We believe that a long lived 
transient is likely to set in. 

It would also be interesting to reexamine our results 
with a dynamical equation allowing for the coexistence 
of loops and infinite strings. We reserved to a future 
publication the presentation of numerical solutions p9| . 
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